This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.
#install.packages("corrplot")
#install.packages("Rfast")
#install.packages("phyloseq")
#install.packages("glmnet")
#install.packages("randcorr")
#install.packages("hrbrthemes")
library(randcorr)
library(ggplot2)
library(hrbrthemes)
library(Rfast)
library(corrplot)
# 1a
create_TC = function(interval, one_time, start_time){
# Create initial matrix
TC = matrix(0, 240, 1)
# Counts for how many ones have been added / not added
on_count = 1
off_count = 0
for(i in 1:240){
# Check if the current time in source is before our start time
if (i <= start_time) {
next
}
# Change value in matrix if it is within the time period to do so
# Otherwise increment the time that one has not been added
# Reset if at the end of a interval
if(on_count <= one_time) {
TC[i] = 1
on_count = on_count + 1
} else{
off_count = off_count + 1
if (off_count == interval - one_time) {
on_count = 1
off_count = 0
}
}
}
return(TC)
}
# 1a
TC1 = create_TC(30, 15, 0)
TC2 = create_TC(45, 20, 20)
TC3 = create_TC(60, 25, 0)
TC4 = create_TC(40, 15, 0)
TC5 = create_TC(40, 20, 0)
TC6 = create_TC(40, 25, 0)
# 1a
standard_TC1 = scale(TC1)
standard_TC2 = scale(TC2)
standard_TC3 = scale(TC3)
standard_TC4 = scale(TC4)
standard_TC5 = scale(TC5)
standard_TC6 = scale(TC6)
full_TC = matrix(c(standard_TC1, standard_TC2, standard_TC3, standard_TC4, standard_TC5, standard_TC6), 240, 6)
#normalised_TC1 = sqrt(sum(TC1^2))
#normalised_TC2 = sqrt(sum(TC2^2))
#normalised_TC3 = sqrt(sum(TC3^2))
#normalised_TC4 = sqrt(sum(TC4^2))
#normalised_TC5 = sqrt(sum(TC5^2))
#normalised_TC6 = sqrt(sum(TC6^2))
#(TC1 / normalised_TC1)[1:50]
#(TC5 / normalised_TC5)[1:50]
#(TC6 / normalised_TC6)[1:50]
# 1a
par(mfrow=c(2,3))
plot(full_TC[1:nrow(full_TC), 1], ylab = "Temporal source 1", pch = 1, type = "line")
plot(full_TC[1:nrow(full_TC), 2], ylab = "Temporal source 2", pch = 1, type = "line")
plot(full_TC[1:nrow(full_TC), 3], ylab = "Temporal source 3", pch = 1, type = "line")
plot(full_TC[1:nrow(full_TC), 4], ylab = "Temporal source 4", pch = 1, type = "line")
plot(full_TC[1:nrow(full_TC), 5], ylab = "Temporal source 5", pch = 1, type = "line")
plot(full_TC[1:nrow(full_TC), 6], ylab = "Temporal source 6", pch = 1, type = "line")

# If we normalise, the values of TC will be same where as standardised values
# Are not the same inter TC eg
# All TC will have 0 values still
# 1b
correlation = cor(full_TC)
corrplot.mixed(correlation, tl.srt = 0)

correlation
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.000000e+00 3.175488e-21 1.690309e-01 0.08606630 -4.536411e-22 0.0860663
[2,] 3.175488e-21 1.000000e+00 -2.857143e-02 0.13093073 1.769200e-20 -0.1309307
[3,] 1.690309e-01 -2.857143e-02 1.000000e+00 0.04364358 4.082770e-21 0.1309307
[4,] 8.606630e-02 1.309307e-01 4.364358e-02 1.00000000 7.745967e-01 0.6000000
[5,] -4.536411e-22 1.769200e-20 4.082770e-21 0.77459667 1.000000e+00 0.7745967
[6,] 8.606630e-02 -1.309307e-01 1.309307e-01 0.60000000 7.745967e-01 1.0000000
# Comment on if on TC correlation or not
# 4/5, 5/6 highly correlated as size of circle is size of correlation and relatively positvely correlated
# 1c
create_tmpSM = function(vert_start,vert_end, horo_start, horo_end) {
tmpSM = matrix(0, 21, 21)
tmpSM[vert_start:vert_end, horo_start:horo_end] = 1
return(tmpSM)
}
# 1c
tmpSM1 = create_tmpSM(2, 6, 2, 6)
tmpSM2 = create_tmpSM(2, 6, 15, 19)
tmpSM3 = create_tmpSM(8, 13, 2, 6)
tmpSM4 = create_tmpSM(8, 13, 15, 19)
tmpSM5 = create_tmpSM(15, 19, 2, 6)
tmpSM6 = create_tmpSM(15, 19, 15, 19)
create_df = function(vert_start,vert_end, horo_start, horo_end) {
df1 = expand.grid(seq_len(21), seq_len(21))
tmp = matrix(0,21,21)
tmp[vert_start:vert_end, horo_start:horo_end] = 1
column = as.vector(tmp)
df1$data = (column)
return(df1)
}
plot1 = create_df(2, 6, 2, 6)
plot2 = create_df(2, 6, 15, 19)
plot3 = create_df(8, 13, 2, 6)
plot4 = create_df(8, 13, 15, 19)
plot5 = create_df(15, 19, 2, 6)
plot6 = create_df(15, 19, 15, 19)
ggplot(plot1, aes(Var1, Var2, fill= data)) +
geom_tile() + scale_fill_gradient(low="blue", high="red") +
theme_ipsum()

ggplot(plot2, aes(Var1, Var2, fill= data)) +
geom_tile() + scale_fill_gradient(low="blue", high="red") +
theme_ipsum()

ggplot(plot3, aes(Var1, Var2, fill= data)) +
geom_tile() + scale_fill_gradient(low="blue", high="red") +
theme_ipsum()

ggplot(plot4, aes(Var1, Var2, fill= data)) +
geom_tile() + scale_fill_gradient(low="blue", high="red") +
theme_ipsum()

ggplot(plot5, aes(Var1, Var2, fill= data)) +
geom_tile() + scale_fill_gradient(low="blue", high="red") +
theme_ipsum()

ggplot(plot6, aes(Var1, Var2, fill= data)) +
geom_tile() + scale_fill_gradient(low="blue", high="red") +
theme_ipsum()

# 1c
tmpSM1 = matrix(tmpSM1, 1, 441)
tmpSM2 = matrix(tmpSM2, 1, 441)
tmpSM3 = matrix(tmpSM3, 1, 441)
tmpSM4 = matrix(tmpSM4, 1, 441)
tmpSM5 = matrix(tmpSM5, 1, 441)
tmpSM6 = matrix(tmpSM6, 1, 441)
tmpSM = matrix(c(tmpSM1, tmpSM2, tmpSM3, tmpSM4, tmpSM5, tmpSM6), 441, 6)
tmpSM = t(tmpSM)
correlation2 = cor(t(tmpSM))
corrplot.mixed(correlation2, upper = "circle")

correlation2
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.00000000 -0.06009615 -0.06623127 -0.06623127 -0.06009615 -0.06009615
[2,] -0.06009615 1.00000000 -0.06623127 -0.06623127 -0.06009615 -0.06009615
[3,] -0.06623127 -0.06623127 1.00000000 -0.07299270 -0.06623127 -0.06623127
[4,] -0.06623127 -0.06623127 -0.07299270 1.00000000 -0.06623127 -0.06623127
[5,] -0.06009615 -0.06009615 -0.06623127 -0.06623127 1.00000000 -0.06009615
[6,] -0.06009615 -0.06009615 -0.06623127 -0.06623127 -0.06009615 1.00000000
# Are theres SM all independent? - yes as not correlated
# Why standardisation not important - If we are looking to replicate it then
# We still just have two main areas - zero and non zero so the scale of non 0
# Doesn't matter as much
#1d
set.seed(30034)
noise_t = rnorm(6 * 240, mean = 0, sd = sqrt(0.25))
noise_s = rnorm(6 * 441, mean = 0, sd = sqrt(0.15))
noise_s = matrix(noise_s, 6, 441)
noise_t = matrix(noise_t, 240, 6)
# 1d
correlation3 = cor(noise_t)
corrplot(correlation3, tl.srt = 0)

correlation4 = cor(t(noise_s))
corrplot(correlation4, tl.srt = 0)

# Check if these two above are correlated
# Highest is like abs(0.12,0.15) so not really that strongly correlated
# 1d
# Both look pretty normally distributed with over 95% of the data covered within -1.96, 1.96
hist(noise_s)

hist(noise_t)

# Both are equal to 1 even with formula variance = 1.96sd, noise_s = 0.945 so pre close, noise_t slight less 0.936 makes sense with graph
#length(noise_s[(0 - 2*sqrt(0.15) <= noise_s) & (noise_s <= 0 + 2*sqrt(0.15))]) / length(noise_s)
#length(noise_t[(0 - 2*sqrt(0.25) <= noise_t) & (noise_t <= 0 + 2*sqrt(0.25))]) / length(noise_t)
# noise s looks way more normally distributed, while t looks kinda skewed.
# do conf int to normal dist
# Check whether fills var of 1.96*sd (0.95 data)
correlation5 = cor(noise_t %*% noise_s)
#corrplot(correlation5, tl.srt = 0)
subset = correlation5[1:10, 1:10]
corrplot(subset, tl.srt = 0)

subset
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1.00000000 -0.08401539 0.08121198 -0.48582752 -0.2122397 -0.58300192 0.7469938 0.1747322 -0.4616947 0.69749810
[2,] -0.08401539 1.00000000 -0.36170166 -0.08651628 0.3311008 -0.65485963 0.1186072 -0.6562960 0.2556765 -0.27880412
[3,] 0.08121198 -0.36170166 1.00000000 -0.61666983 -0.2739593 0.05563854 0.4470863 0.6308545 -0.6288561 -0.05763502
[4,] -0.48582752 -0.08651628 -0.61666983 1.00000000 0.2248447 0.47069260 -0.5268018 -0.4475564 0.4159525 -0.37703467
[5,] -0.21223971 0.33110081 -0.27395932 0.22484470 1.0000000 0.24773797 -0.1582580 -0.5139208 0.6008665 -0.70374353
[6,] -0.58300192 -0.65485963 0.05563854 0.47069260 0.2477380 1.00000000 -0.6314982 0.2674007 0.3960874 -0.40361348
[7,] 0.74699385 0.11860722 0.44708627 -0.52680181 -0.1582580 -0.63149817 1.0000000 0.2406842 -0.5654764 0.25340041
[8,] 0.17473222 -0.65629601 0.63085453 -0.44755642 -0.5139208 0.26740068 0.2406842 1.0000000 -0.2328050 0.37238708
[9,] -0.46169469 0.25567646 -0.62885608 0.41595249 0.6008665 0.39608744 -0.5654764 -0.2328050 1.0000000 -0.39341424
[10,] 0.69749810 -0.27880412 -0.05763502 -0.37703467 -0.7037435 -0.40361348 0.2534004 0.3723871 -0.3934142 1.00000000
#corrplot(correlation5, tl.srt = 0)[1:10, 1:10]
# Definitely some correlation eg higher covariance values such as 0.7
#1e
#Can full_tc%*% noise_S and vice versa exist and what happens to them as if we keep we cannot fit model?
X = (full_TC + noise_t) %*% (tmpSM + noise_s)
scaled_X = matrix(NA, 240, 441)
for (i in 1:441) {
scaled_X[,i] = scale(X[,i])
}
# Plot graph of sample Time series
set.seed(1)
mysample <- sample(nrow(X), 100)
mysample
[1] 68 167 129 162 215 43 14 210 187 51 225 85 21 106 182 74 7 73 79 213 37 105 217 110 165 34 227 126 89 172 207
[32] 33 84 163 70 230 42 166 111 148 156 20 44 121 87 206 197 40 211 25 119 122 39 179 204 134 24 160 234 130 45 146
[63] 22 115 104 195 175 176 103 75 13 193 212 23 218 208 192 29 141 150 200 108 48 209 180 149 31 102 145 223 196 83 118
[94] 90 158 107 64 164 60 231
sample = X[1:nrow(X), mysample]
#for (j in 1:10) {
# if (j == 1) {
# plot(sample[1:nrow(X), (i - 1)*10 + j], type = "line")
# }
# else {
# if (j == 2) {
# lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'red')
# }
#
# if (j == 3) {
# lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'blue')
# }
#
# if (j == 4) {
# lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'green')
# }
#
# if (j == 5) {
# lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'yellow')
# }
# if (j == 6) {
# lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'darkorange')
# }
#
# if (j == 7) {
# lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'cyan')
# }
# if (j == 8) {
# lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'brown')
# }
# if (j == 9) {
# lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'grey')
# }
#
# if (j == 10) {
# lines(sample[1:nrow(X), (i - 1)*10 + j], col = 'beige')
# }
#
# }
#
#}
# What info does the plot of vars give
variances = colVars(X)
plot(variances, type = "line")

NA
NA
NA
# 2a
A_LSR = solve(t(full_TC) %*% full_TC) %*% t(full_TC) %*% scaled_X
D_LSR = scaled_X%*% abs(t(A_LSR))
# Plot and heatmap together
plot(D_LSR[1:nrow(D_LSR),1], type = "line")

# produce legend for heatmap
#heatmap(matrix(A[1,1:ncol(A_LSR)], 21, 21), ncol, Colv = NA, Rowv = NA, scale="column")
for (num in 1:6) {
df = as.data.frame(matrix(rep(0, 441 * 3), 441, 3))
count = 1
col = abs(A_LSR[num,])
for (i in 1:21) {
for (j in 1:21) {
df[count, 1] = i
df[count, 2] = j
df[count, 3] = col[(i-1) * 21 + j]
count = count + 1
}
}
ggplot(df, aes(V1, V2, fill= V3)) +
geom_tile() + scale_fill_gradient(low="blue", high="red") +
theme_ipsum()
}
# WHy linear relo bet 3rd col d_lsr and 30th of X
# Secondly if you look at the slice number 3 build up by third SM. Now read this 3rd SM column-wise, you will find that 30th (2nd column 9th row =30th pixel) pixel position was filled by this 3rd SM and therefore third TC is the only time course that constructs 30th column of X and this is the TC that you retrieve in 3rd Dlsr. There might be some variance in the scatter plot but you will be able to capture visually the linear trend between 3rd Dlsr and 30th X.
#linear relo -
plot(D_LSR[1:nrow(D_LSR), 3], scaled_X[1:nrow(scaled_X), 30], xlab = "3rd column of D(LSR)", ylab = "30th column of scaled X")

# No linear relo
plot(D_LSR[1:nrow(D_LSR), 4], scaled_X[1:nrow(scaled_X), 30], xlab = "3rd column of D(LSR)", ylab = "30th column of scaled X")

# Gets the sum of the maximum absolute values in the correlation matrix
max_abs_corr_sum = function(mat1, mat2) {
corre_mat = cor(mat1, mat2)
c = get_max_correlations(corre_mat)
return(sum(c))
}
# Returns the values of the maximum absolute values for each row
get_max_correlations = function(cor_mat) {
c_var = c(rep(0,6))
abs_mat = abs(cor_mat)
cols = max.col(abs_mat)
for (i in 1:6) {
c_var[i] = cor_mat[i, cols[i]]
}
return(c_var)
}
#2b
V = 441
MSE = matrix(0, 1, 2)
for (num in seq(from = 0, to = 1, by = 0.01)) {
lamda = num
lamda_hat = lamda * V
A_RR = solve(t(full_TC) %*% full_TC + diag(nrow(t(full_TC))) * lamda_hat) %*% t(full_TC) %*% scaled_X
mse = (sum((scaled_X - full_TC %*% abs(A_RR))^2)) + lamda_hat * sum(A_RR^2)
MSE = rbind(MSE, c(lamda, mse))
}
plot(MSE[2:nrow(MSE), 1], MSE[2:nrow(MSE), 2])

# From the graph it looks like best value of lamda is 1
A_RR = solve(t(full_TC) %*% full_TC + diag(nrow(t(full_TC))) * 1* V) %*% t(full_TC) %*% scaled_X
D_RR = scaled_X%*% abs(t(A_RR))
c_tlsr = max_abs_corr_sum(full_TC, D_LSR)
c_trr = max_abs_corr_sum(full_TC, D_RR)
# 4.407241
c_tlsr
[1] 4.506323
# 4.47653
c_trr
[1] 4.682376
c_trr - c_tlsr
[1] 0.1760527
lamda = 1000
lamda_hat = lamda * V
A_1000 = solve(t(full_TC) %*% full_TC + diag(nrow(t(full_TC))) * lamda_hat) %*% t(full_TC) %*% scaled_X
# As A_LSR decrease, A_RR values decrease
plot(abs(A_1000[1:nrow(A_1000), 1]), abs(A_LSR[1:nrow(A_LSR), 1]), xlab = "A_RR", ylab = "A_LSR")

create_noise_matrix = function(seed) {
set.seed(seed = seed)
noise_t = rnorm(6 * 240, mean = 0, sd = sqrt(0.25))
noise_s = rnorm(6 * 441, mean = 0, sd = sqrt(0.15))
noise_s = matrix(noise_s, 6, 441)
noise_t = matrix(noise_t, 240, 6)
X = (full_TC + noise_t) %*% (tmpSM + noise_s)
new_X = scale(X)
return(new_X)
}
noise_1 = create_noise_matrix(30035)
noise_2 = create_noise_matrix(30037)
noise_3 = create_noise_matrix(30038)
noise_4 = create_noise_matrix(30036)
noise_5 = create_noise_matrix(30039)
noise_6 = create_noise_matrix(30040)
noise_7 = create_noise_matrix(30041)
noise_8 = create_noise_matrix(30042)
noise_9 = create_noise_matrix(30043)
noise_10 = create_noise_matrix(30044)
noise_df = list(noise_1, noise_9, noise_8, noise_7, noise_6, noise_5,
noise_2, noise_3, noise_4, noise_10
)
#for(mat in noise_df) {
# print(mat[1:5, 1:3]
#}
#2c
LR_MSE = matrix(NA, 21, 10)
N = 240
x1 = 21
x2 = 21
nsrcs = 6
V = x1 * x2
i = 1
j = 1
step <- 1/(norm(full_TC %*% t(full_TC)) * 1.1)
for (num in seq(0,1,0.05)) {
for (x_rep in noise_df) {
thr <- num*N*step
Ao <- matrix(0, nsrcs, 1)
A <- matrix(0, nsrcs, 1)
Alr <- matrix(0, nsrcs, x1*x2)
for (k in 1:(x1*x2)) {
A <- Ao+step*(t(full_TC) %*% (x_rep[,k]-(full_TC%*%Ao)))
A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
for (rep in 1:10) {
Ao <- A
A <- Ao+step * (t(full_TC)%*%(x_rep[,k]-(full_TC%*%Ao)))
A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
}
Alr[,k] <- A
}
Dlr = x_rep %*% abs(t(Alr))
LR_MSE[i,j] = sum(sum((x_rep-Dlr%*%abs(Alr))^2)) / (N*V)
j = j + 1
}
j = 1
i = i + 1
}
LR_avgs = rowMeans(LR_MSE)
plot( seq(0,1, 0.05), LR_avgs, ylab = "MSE(LR)")

# Min value is 0.580 when p = 0.55, at same value it diverges
# Since first index is at 0, we want the 12th value between 0 and 1 which is 0.6
# So we need to index the 13th value
index = which.min(LR_avgs)
rho = (index - 1) * 0.05
# 2d
thr <- rho*N*step
Ao <- matrix(0, nsrcs, 1)
A <- matrix(0, nsrcs, 1)
Alr <- matrix(0, nsrcs, x1*x2)
for (k in 1:(x1*x2)) {
A <- Ao+step*(t(full_TC) %*% (scaled_X[,k]-(full_TC%*%Ao)))
A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
for (rep in 1:10) {
Ao <- A
A <- Ao+step * (t(full_TC)%*%(scaled_X[,k]-(full_TC%*%Ao)))
A <- (1/(1+thr)) * (sign(A)*pmax(replicate(nsrcs, 0), abs(A)-thr))
}
Alr[,k] <- A
}
Dlr = scaled_X %*% abs(t(Alr))
ctrr = max_abs_corr_sum(full_TC, D_RR)
csrr = max_abs_corr_sum(t(tmpSM), t(A_RR))
ctlr = max_abs_corr_sum(full_TC, Dlr)
cslr= max_abs_corr_sum(t(tmpSM), t(Alr))
# Both are true
ctlr - ctrr #= 0.7413487
[1] 0.6743216
cslr - csrr #= 0.6595151
[1] 0.6778317
ctlr #= 5.217879
[1] 5.356698
ctrr #= 4.47653
[1] 4.682376
cslr #= 2.953101
[1] 2.937037
csrr #= 2.293586
[1] 2.259206
plot(D_RR[, 1:4])

plot(Dlr[, 1:4])

plot(Alr[1:4,])

plot(A_RR[1:4, ])
par(mfrow = c(1,2))

for (num in 1:6) {
plot(D_RR[, 6], type = "line")
plot(Dlr[, 6], type = "line")
}

plot type 'line' will be truncated to first character

plot type 'line' will be truncated to first character

plot type 'line' will be truncated to first character

plot type 'line' will be truncated to first character

plot type 'line' will be truncated to first character

df = as.data.frame(matrix(rep(0, 441 * 3), 441, 3))
count = 1
vec = abs(A_RR[num,])
for (num in 1:6) {
for (i in 1:21) {
for (j in 1:21) {
df[count, 1] = i
df[count, 2] = j
df[count, 3] = vec[(i-1) * 21 + j]
count = count + 1
}
}
}
df = as.data.frame(matrix(rep(0, 441 * 3), 441, 3))
count = 1
vec = abs(Alr[num,])
for (num in 1:6) {
for (i in 1:21) {
for (j in 1:21) {
df[count, 1] = i
df[count, 2] = j
df[count, 3] = vec[(i-1) * 21 + j]
count = count + 1
}
}
}
ggplot(df, aes(V1, V2, fill= V3)) +
geom_tile() + scale_fill_gradient(low="blue", high="red") +
theme_ipsum()

par(mfrow=c(1,2))
# 2e
#svd(full_TC)
plot(svd(full_TC)$d, ylab = "Eigenvalue")
#PC 6 smallest
plot(svd(full_TC)$u[1:240, 6])

plot(full_TC[1:240, 6])
plot(svd(full_TC)$u[1:240, 5])

plot(full_TC[1:240, 5])
plot(svd(full_TC)$u[1:240, 4])

plot(full_TC[1:240, 4])
plot(svd(full_TC)$u[1:240, 3])

plot(full_TC[1:240, 3])
plot(svd(full_TC)$u[1:240, 2])

plot(full_TC[1:240, 2])
plot(svd(full_TC)$u[1:240, 1])

plot(full_TC[1:240, 1])

NA
NA
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKVGhpcyBpcyBhbiBbUiBNYXJrZG93bl0oaHR0cDovL3JtYXJrZG93bi5yc3R1ZGlvLmNvbSkgTm90ZWJvb2suIFdoZW4geW91IGV4ZWN1dGUgY29kZSB3aXRoaW4gdGhlIG5vdGVib29rLCB0aGUgcmVzdWx0cyBhcHBlYXIgYmVuZWF0aCB0aGUgY29kZS4gCgpUcnkgZXhlY3V0aW5nIHRoaXMgY2h1bmsgYnkgY2xpY2tpbmcgdGhlICpSdW4qIGJ1dHRvbiB3aXRoaW4gdGhlIGNodW5rIG9yIGJ5IHBsYWNpbmcgeW91ciBjdXJzb3IgaW5zaWRlIGl0IGFuZCBwcmVzc2luZyAqQ21kK1NoaWZ0K0VudGVyKi4gCmBgYHtyfQojaW5zdGFsbC5wYWNrYWdlcygiY29ycnBsb3QiKQojaW5zdGFsbC5wYWNrYWdlcygiUmZhc3QiKQojaW5zdGFsbC5wYWNrYWdlcygicGh5bG9zZXEiKQojaW5zdGFsbC5wYWNrYWdlcygiZ2xtbmV0IikKI2luc3RhbGwucGFja2FnZXMoInJhbmRjb3JyIikKI2luc3RhbGwucGFja2FnZXMoImhyYnJ0aGVtZXMiKQpsaWJyYXJ5KHJhbmRjb3JyKQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkoaHJicnRoZW1lcykKbGlicmFyeShSZmFzdCkKbGlicmFyeShjb3JycGxvdCkKYGBgCgoKYGBge3J9CiMgMWEKY3JlYXRlX1RDID0gZnVuY3Rpb24oaW50ZXJ2YWwsIG9uZV90aW1lLCBzdGFydF90aW1lKXsKICAKICAjIENyZWF0ZSBpbml0aWFsIG1hdHJpeAogIFRDID0gbWF0cml4KDAsIDI0MCwgMSkKICAKICAjIENvdW50cyBmb3IgaG93IG1hbnkgb25lcyBoYXZlIGJlZW4gYWRkZWQgLyBub3QgYWRkZWQKICBvbl9jb3VudCA9IDEKICBvZmZfY291bnQgPSAwCiAgCiAgZm9yKGkgaW4gMToyNDApewogICAgCiAgICAjIENoZWNrIGlmIHRoZSBjdXJyZW50IHRpbWUgaW4gc291cmNlIGlzIGJlZm9yZSBvdXIgc3RhcnQgdGltZQogICAgaWYgKGkgPD0gc3RhcnRfdGltZSkgewogICAgICBuZXh0CiAgICB9CiAgICAKICAgICMgQ2hhbmdlIHZhbHVlIGluIG1hdHJpeCBpZiBpdCBpcyB3aXRoaW4gdGhlIHRpbWUgcGVyaW9kIHRvIGRvIHNvIAogICAgIyBPdGhlcndpc2UgaW5jcmVtZW50IHRoZSB0aW1lIHRoYXQgb25lIGhhcyBub3QgYmVlbiBhZGRlZAogICAgIyBSZXNldCBpZiBhdCB0aGUgZW5kIG9mIGEgaW50ZXJ2YWwKICAgIGlmKG9uX2NvdW50IDw9IG9uZV90aW1lKSB7CiAgICAgIFRDW2ldID0gMQogICAgICBvbl9jb3VudCA9IG9uX2NvdW50ICsgMQogICAgfSBlbHNlewogICAgICBvZmZfY291bnQgPSBvZmZfY291bnQgKyAxCiAgICAgIGlmIChvZmZfY291bnQgPT0gaW50ZXJ2YWwgLSBvbmVfdGltZSkgewogICAgICAgIG9uX2NvdW50ID0gMQogICAgICAgIG9mZl9jb3VudCA9IDAKICAgICAgfQogIAogICAgfQogIAogIH0KICByZXR1cm4oVEMpCn0KYGBgCgoKYGBge3J9CiMgMWEKVEMxID0gY3JlYXRlX1RDKDMwLCAxNSwgMCkKVEMyID0gY3JlYXRlX1RDKDQ1LCAyMCwgMjApClRDMyA9IGNyZWF0ZV9UQyg2MCwgMjUsIDApClRDNCA9IGNyZWF0ZV9UQyg0MCwgMTUsIDApClRDNSA9IGNyZWF0ZV9UQyg0MCwgMjAsIDApClRDNiA9IGNyZWF0ZV9UQyg0MCwgMjUsIDApCgpgYGAKCmBgYHtyfQojIDFhCnN0YW5kYXJkX1RDMSA9IHNjYWxlKFRDMSkKc3RhbmRhcmRfVEMyID0gc2NhbGUoVEMyKQpzdGFuZGFyZF9UQzMgPSBzY2FsZShUQzMpCnN0YW5kYXJkX1RDNCA9IHNjYWxlKFRDNCkKc3RhbmRhcmRfVEM1ID0gc2NhbGUoVEM1KQpzdGFuZGFyZF9UQzYgPSBzY2FsZShUQzYpCmZ1bGxfVEMgPSBtYXRyaXgoYyhzdGFuZGFyZF9UQzEsIHN0YW5kYXJkX1RDMiwgc3RhbmRhcmRfVEMzLCBzdGFuZGFyZF9UQzQsIHN0YW5kYXJkX1RDNSwgc3RhbmRhcmRfVEM2KSwgMjQwLCA2KQoKI25vcm1hbGlzZWRfVEMxID0gc3FydChzdW0oVEMxXjIpKQojbm9ybWFsaXNlZF9UQzIgPSBzcXJ0KHN1bShUQzJeMikpCiNub3JtYWxpc2VkX1RDMyA9IHNxcnQoc3VtKFRDM14yKSkKI25vcm1hbGlzZWRfVEM0ID0gc3FydChzdW0oVEM0XjIpKQojbm9ybWFsaXNlZF9UQzUgPSBzcXJ0KHN1bShUQzVeMikpCiNub3JtYWxpc2VkX1RDNiA9IHNxcnQoc3VtKFRDNl4yKSkKCiMoVEMxIC8gbm9ybWFsaXNlZF9UQzEpWzE6NTBdCiMoVEM1IC8gbm9ybWFsaXNlZF9UQzUpWzE6NTBdCiMoVEM2IC8gbm9ybWFsaXNlZF9UQzYpWzE6NTBdCmBgYAoKCmBgYHtyfQojIDFhCnBhcihtZnJvdz1jKDIsMykpCgpwbG90KGZ1bGxfVENbMTpucm93KGZ1bGxfVEMpLCAxXSwgeWxhYiA9ICJUZW1wb3JhbCBzb3VyY2UgMSIsIHBjaCA9IDEsIHR5cGUgPSAibGluZSIpCnBsb3QoZnVsbF9UQ1sxOm5yb3coZnVsbF9UQyksIDJdLCB5bGFiID0gIlRlbXBvcmFsIHNvdXJjZSAyIiwgcGNoID0gMSwgdHlwZSA9ICJsaW5lIikKcGxvdChmdWxsX1RDWzE6bnJvdyhmdWxsX1RDKSwgM10sIHlsYWIgPSAiVGVtcG9yYWwgc291cmNlIDMiLCBwY2ggPSAxLCB0eXBlID0gImxpbmUiKQpwbG90KGZ1bGxfVENbMTpucm93KGZ1bGxfVEMpLCA0XSwgeWxhYiA9ICJUZW1wb3JhbCBzb3VyY2UgNCIsIHBjaCA9IDEsIHR5cGUgPSAibGluZSIpCnBsb3QoZnVsbF9UQ1sxOm5yb3coZnVsbF9UQyksIDVdLCB5bGFiID0gIlRlbXBvcmFsIHNvdXJjZSA1IiwgcGNoID0gMSwgdHlwZSA9ICJsaW5lIikKcGxvdChmdWxsX1RDWzE6bnJvdyhmdWxsX1RDKSwgNl0sIHlsYWIgPSAiVGVtcG9yYWwgc291cmNlIDYiLCBwY2ggPSAxLCB0eXBlID0gImxpbmUiKQoKCgojIElmIHdlIG5vcm1hbGlzZSwgdGhlIHZhbHVlcyBvZiBUQyB3aWxsIGJlIHNhbWUgd2hlcmUgYXMgc3RhbmRhcmRpc2VkIHZhbHVlcyAKIyBBcmUgbm90IHRoZSBzYW1lIGludGVyIFRDIGVnCiMgQWxsIFRDIHdpbGwgaGF2ZSAwIHZhbHVlcyBzdGlsbApgYGAKCmBgYHtyfQojIDFiCgpjb3JyZWxhdGlvbiA9IGNvcihmdWxsX1RDKQpjb3JycGxvdC5taXhlZChjb3JyZWxhdGlvbiwgdGwuc3J0ID0gMCkKCmNvcnJlbGF0aW9uCiMgQ29tbWVudCBvbiBpZiBvbiBUQyBjb3JyZWxhdGlvbiBvciBub3QKIyA0LzUsIDUvNiBoaWdobHkgY29ycmVsYXRlZCBhcyBzaXplIG9mIGNpcmNsZSBpcyBzaXplIG9mIGNvcnJlbGF0aW9uIGFuZCByZWxhdGl2ZWx5IHBvc2l0dmVseSBjb3JyZWxhdGVkCmBgYAoKYGBge3J9CiMgMWMKY3JlYXRlX3RtcFNNID0gZnVuY3Rpb24odmVydF9zdGFydCx2ZXJ0X2VuZCwgaG9yb19zdGFydCwgaG9yb19lbmQpIHsKICB0bXBTTSA9IG1hdHJpeCgwLCAyMSwgMjEpCiAgdG1wU01bdmVydF9zdGFydDp2ZXJ0X2VuZCwgaG9yb19zdGFydDpob3JvX2VuZF0gPSAxCiAgcmV0dXJuKHRtcFNNKQp9CmBgYAoKCmBgYHtyfQojIDFjCnRtcFNNMSA9IGNyZWF0ZV90bXBTTSgyLCA2LCAyLCA2KQp0bXBTTTIgPSBjcmVhdGVfdG1wU00oMiwgNiwgMTUsIDE5KQp0bXBTTTMgPSBjcmVhdGVfdG1wU00oOCwgMTMsIDIsIDYpCnRtcFNNNCA9IGNyZWF0ZV90bXBTTSg4LCAxMywgMTUsIDE5KQp0bXBTTTUgPSBjcmVhdGVfdG1wU00oMTUsIDE5LCAyLCA2KQp0bXBTTTYgPSBjcmVhdGVfdG1wU00oMTUsIDE5LCAxNSwgMTkpCmBgYAoKYGBge3J9CgpjcmVhdGVfZGYgPSBmdW5jdGlvbih2ZXJ0X3N0YXJ0LHZlcnRfZW5kLCBob3JvX3N0YXJ0LCBob3JvX2VuZCkgewogIGRmMSA9IGV4cGFuZC5ncmlkKHNlcV9sZW4oMjEpLCBzZXFfbGVuKDIxKSkKICB0bXAgPSBtYXRyaXgoMCwyMSwyMSkKICB0bXBbdmVydF9zdGFydDp2ZXJ0X2VuZCwgaG9yb19zdGFydDpob3JvX2VuZF0gPSAxCiAgY29sdW1uID0gYXMudmVjdG9yKHRtcCkKICBkZjEkZGF0YSA9IChjb2x1bW4pCiAgcmV0dXJuKGRmMSkKfQoKcGxvdDEgPSBjcmVhdGVfZGYoMiwgNiwgMiwgNikKcGxvdDIgPSBjcmVhdGVfZGYoMiwgNiwgMTUsIDE5KQpwbG90MyA9IGNyZWF0ZV9kZig4LCAxMywgMiwgNikKcGxvdDQgPSBjcmVhdGVfZGYoOCwgMTMsIDE1LCAxOSkKcGxvdDUgPSBjcmVhdGVfZGYoMTUsIDE5LCAyLCA2KQpwbG90NiA9IGNyZWF0ZV9kZigxNSwgMTksIDE1LCAxOSkKCmdncGxvdChwbG90MSwgYWVzKFZhcjEsIFZhcjIsIGZpbGw9IGRhdGEpKSArIAogIGdlb21fdGlsZSgpICsgc2NhbGVfZmlsbF9ncmFkaWVudChsb3c9ImJsdWUiLCBoaWdoPSJyZWQiKSArCiAgdGhlbWVfaXBzdW0oKQoKZ2dwbG90KHBsb3QyLCBhZXMoVmFyMSwgVmFyMiwgZmlsbD0gZGF0YSkpICsgCiAgZ2VvbV90aWxlKCkgKyBzY2FsZV9maWxsX2dyYWRpZW50KGxvdz0iYmx1ZSIsIGhpZ2g9InJlZCIpICsKICB0aGVtZV9pcHN1bSgpCgpnZ3Bsb3QocGxvdDMsIGFlcyhWYXIxLCBWYXIyLCBmaWxsPSBkYXRhKSkgKyAKICBnZW9tX3RpbGUoKSArIHNjYWxlX2ZpbGxfZ3JhZGllbnQobG93PSJibHVlIiwgaGlnaD0icmVkIikgKwogIHRoZW1lX2lwc3VtKCkKCmdncGxvdChwbG90NCwgYWVzKFZhcjEsIFZhcjIsIGZpbGw9IGRhdGEpKSArIAogIGdlb21fdGlsZSgpICsgc2NhbGVfZmlsbF9ncmFkaWVudChsb3c9ImJsdWUiLCBoaWdoPSJyZWQiKSArCiAgdGhlbWVfaXBzdW0oKQoKZ2dwbG90KHBsb3Q1LCBhZXMoVmFyMSwgVmFyMiwgZmlsbD0gZGF0YSkpICsgCiAgZ2VvbV90aWxlKCkgKyBzY2FsZV9maWxsX2dyYWRpZW50KGxvdz0iYmx1ZSIsIGhpZ2g9InJlZCIpICsKICB0aGVtZV9pcHN1bSgpCgpnZ3Bsb3QocGxvdDYsIGFlcyhWYXIxLCBWYXIyLCBmaWxsPSBkYXRhKSkgKyAKICBnZW9tX3RpbGUoKSArIHNjYWxlX2ZpbGxfZ3JhZGllbnQobG93PSJibHVlIiwgaGlnaD0icmVkIikgKwogIHRoZW1lX2lwc3VtKCkKCmBgYApgYGB7cn0KIyAxYwp0bXBTTTEgPSBtYXRyaXgodG1wU00xLCAxLCA0NDEpCnRtcFNNMiA9IG1hdHJpeCh0bXBTTTIsIDEsIDQ0MSkKdG1wU00zID0gbWF0cml4KHRtcFNNMywgMSwgNDQxKQp0bXBTTTQgPSBtYXRyaXgodG1wU000LCAxLCA0NDEpCnRtcFNNNSA9IG1hdHJpeCh0bXBTTTUsIDEsIDQ0MSkKdG1wU002ID0gbWF0cml4KHRtcFNNNiwgMSwgNDQxKQp0bXBTTSA9IG1hdHJpeChjKHRtcFNNMSwgdG1wU00yLCB0bXBTTTMsIHRtcFNNNCwgdG1wU001LCB0bXBTTTYpLCA0NDEsIDYpCnRtcFNNID0gdCh0bXBTTSkKCmNvcnJlbGF0aW9uMiA9IGNvcih0KHRtcFNNKSkKY29ycnBsb3QubWl4ZWQoY29ycmVsYXRpb24yLCB1cHBlciA9ICJjaXJjbGUiKQpjb3JyZWxhdGlvbjIKIyBBcmUgdGhlcmVzIFNNIGFsbCBpbmRlcGVuZGVudD8gLSB5ZXMgYXMgbm90IGNvcnJlbGF0ZWQKIyBXaHkgc3RhbmRhcmRpc2F0aW9uIG5vdCBpbXBvcnRhbnQgLSBJZiB3ZSBhcmUgbG9va2luZyB0byByZXBsaWNhdGUgaXQgdGhlbiAKIyBXZSBzdGlsbCBqdXN0IGhhdmUgdHdvIG1haW4gYXJlYXMgLSB6ZXJvIGFuZCBub24gemVybyBzbyB0aGUgc2NhbGUgb2Ygbm9uIDAKIyBEb2Vzbid0IG1hdHRlciBhcyBtdWNoCmBgYAoKCmBgYHtyfQojMWQKc2V0LnNlZWQoMzAwMzQpCm5vaXNlX3QgPSBybm9ybSg2ICogMjQwLCBtZWFuID0gMCwgc2QgPSBzcXJ0KDAuMjUpKQpub2lzZV9zID0gcm5vcm0oNiAqIDQ0MSwgbWVhbiA9IDAsIHNkID0gc3FydCgwLjE1KSkKCgpub2lzZV9zID0gbWF0cml4KG5vaXNlX3MsIDYsIDQ0MSkKCm5vaXNlX3QgPSBtYXRyaXgobm9pc2VfdCwgMjQwLCA2KQpgYGAKCgpgYGB7cn0KIyAxZApjb3JyZWxhdGlvbjMgPSBjb3Iobm9pc2VfdCkKY29ycnBsb3QoY29ycmVsYXRpb24zLCB0bC5zcnQgPSAwKQoKY29ycmVsYXRpb240ID0gY29yKHQobm9pc2VfcykpCmNvcnJwbG90KGNvcnJlbGF0aW9uNCwgdGwuc3J0ID0gMCkKIyBDaGVjayBpZiB0aGVzZSB0d28gYWJvdmUgYXJlIGNvcnJlbGF0ZWQKIyBIaWdoZXN0IGlzIGxpa2UgYWJzKDAuMTIsMC4xNSkgc28gbm90IHJlYWxseSB0aGF0IHN0cm9uZ2x5IGNvcnJlbGF0ZWQKYGBgCgoKYGBge3J9CiMgMWQKIyBCb3RoIGxvb2sgcHJldHR5IG5vcm1hbGx5IGRpc3RyaWJ1dGVkIHdpdGggb3ZlciA5NSUgb2YgdGhlIGRhdGEgY292ZXJlZCB3aXRoaW4gLTEuOTYsIDEuOTYKaGlzdChub2lzZV9zKQpoaXN0KG5vaXNlX3QpCgogICMgQm90aCBhcmUgZXF1YWwgdG8gMSBldmVuIHdpdGggZm9ybXVsYSB2YXJpYW5jZSA9IDEuOTZzZCwgbm9pc2VfcyA9IDAuOTQ1IHNvIHByZSBjbG9zZSwgbm9pc2VfdCBzbGlnaHQgbGVzcyAwLjkzNiBtYWtlcyBzZW5zZSB3aXRoIGdyYXBoCiNsZW5ndGgobm9pc2Vfc1soMCAtIDIqc3FydCgwLjE1KSA8PSBub2lzZV9zKSAmIChub2lzZV9zIDw9IDAgKyAyKnNxcnQoMC4xNSkpXSkgLyBsZW5ndGgobm9pc2VfcykKI2xlbmd0aChub2lzZV90WygwIC0gMipzcXJ0KDAuMjUpIDw9IG5vaXNlX3QpICYgKG5vaXNlX3QgPD0gMCArIDIqc3FydCgwLjI1KSldKSAvIGxlbmd0aChub2lzZV90KQoKIyBub2lzZSBzIGxvb2tzIHdheSBtb3JlIG5vcm1hbGx5IGRpc3RyaWJ1dGVkLCB3aGlsZSB0IGxvb2tzIGtpbmRhIHNrZXdlZC4KIyBkbyBjb25mIGludCB0byBub3JtYWwgZGlzdAojIENoZWNrIHdoZXRoZXIgZmlsbHMgdmFyIG9mIDEuOTYqc2QgKDAuOTUgZGF0YSkKCmNvcnJlbGF0aW9uNSA9IGNvcihub2lzZV90ICUqJSBub2lzZV9zKQojY29ycnBsb3QoY29ycmVsYXRpb241LCB0bC5zcnQgPSAwKQpzdWJzZXQgPSBjb3JyZWxhdGlvbjVbMToxMCwgMToxMF0KY29ycnBsb3Qoc3Vic2V0LCB0bC5zcnQgPSAwKQpzdWJzZXQKI2NvcnJwbG90KGNvcnJlbGF0aW9uNSwgdGwuc3J0ID0gMClbMToxMCwgMToxMF0KIyBEZWZpbml0ZWx5IHNvbWUgY29ycmVsYXRpb24gZWcgaGlnaGVyIGNvdmFyaWFuY2UgdmFsdWVzIHN1Y2ggYXMgMC43CgpgYGAKCgpgYGB7cn0KIzFlCgojQ2FuIGZ1bGxfdGMlKiUgbm9pc2VfUyBhbmQgdmljZSB2ZXJzYSBleGlzdCBhbmQgd2hhdCBoYXBwZW5zIHRvIHRoZW0gYXMgaWYgd2Uga2VlcCB3ZSBjYW5ub3QgZml0IG1vZGVsPwpYID0gKGZ1bGxfVEMgKyBub2lzZV90KSAlKiUgKHRtcFNNICsgbm9pc2VfcykKc2NhbGVkX1ggPSBtYXRyaXgoTkEsIDI0MCwgNDQxKQpmb3IgKGkgaW4gMTo0NDEpIHsKICBzY2FsZWRfWFssaV0gPSBzY2FsZShYWyxpXSkKfQoKCiMgUGxvdCBncmFwaCBvZiBzYW1wbGUgVGltZSBzZXJpZXMKc2V0LnNlZWQoMSkKbXlzYW1wbGUgPC0gc2FtcGxlKG5yb3coWCksIDEwMCkKbXlzYW1wbGUKc2FtcGxlID0gWFsxOm5yb3coWCksIG15c2FtcGxlXQoKI2ZvciAoaiBpbiAxOjEwKSB7CiMgIGlmIChqID09IDEpIHsKIyAgICBwbG90KHNhbXBsZVsxOm5yb3coWCksIChpIC0gMSkqMTAgKyBqXSwgdHlwZSA9ICJsaW5lIikKIyAgfQojICBlbHNlIHsKIyAgICBpZiAoaiA9PSAyKSB7CiMgICAgICBsaW5lcyhzYW1wbGVbMTpucm93KFgpLCAoaSAtIDEpKjEwICsgal0sIGNvbCA9ICdyZWQnKQojICAgIH0KIyAgICAKIyAgICBpZiAoaiA9PSAzKSB7CiMgICAgICBsaW5lcyhzYW1wbGVbMTpucm93KFgpLCAoaSAtIDEpKjEwICsgal0sIGNvbCA9ICdibHVlJykKIyAgICB9CiMgICAgCiMgICAgaWYgKGogPT0gNCkgewojICAgICAgbGluZXMoc2FtcGxlWzE6bnJvdyhYKSwgKGkgLSAxKSoxMCArIGpdLCBjb2wgPSAnZ3JlZW4nKQojICAgIH0KIyAgICAKIyAgICBpZiAoaiA9PSA1KSB7CiMgICAgICBsaW5lcyhzYW1wbGVbMTpucm93KFgpLCAoaSAtIDEpKjEwICsgal0sIGNvbCA9ICd5ZWxsb3cnKQojICAgIH0KIyAgICBpZiAoaiA9PSA2KSB7CiMgICAgICBsaW5lcyhzYW1wbGVbMTpucm93KFgpLCAoaSAtIDEpKjEwICsgal0sIGNvbCA9ICdkYXJrb3JhbmdlJykKIyAgICB9CiMgICAgCiMgICAgaWYgKGogPT0gNykgewojICAgICAgbGluZXMoc2FtcGxlWzE6bnJvdyhYKSwgKGkgLSAxKSoxMCArIGpdLCBjb2wgPSAnY3lhbicpCiMgICAgfQojICAgIGlmIChqID09IDgpIHsKIyAgICAgIGxpbmVzKHNhbXBsZVsxOm5yb3coWCksIChpIC0gMSkqMTAgKyBqXSwgY29sID0gJ2Jyb3duJykKIyAgICB9CiMgICAgaWYgKGogPT0gOSkgewojICAgICAgbGluZXMoc2FtcGxlWzE6bnJvdyhYKSwgKGkgLSAxKSoxMCArIGpdLCBjb2wgPSAnZ3JleScpCiMgICAgfQojICAgIAojICAgIGlmIChqID09IDEwKSB7CiMgICAgICBsaW5lcyhzYW1wbGVbMTpucm93KFgpLCAoaSAtIDEpKjEwICsgal0sIGNvbCA9ICdiZWlnZScpCiMgICAgfQojICAgIAojICB9CiMgIAojfQoKIyBXaGF0IGluZm8gZG9lcyB0aGUgcGxvdCBvZiB2YXJzIGdpdmUKdmFyaWFuY2VzID0gY29sVmFycyhYKQpwbG90KHZhcmlhbmNlcywgdHlwZSA9ICJsaW5lIikKCgoKYGBgCgoKYGBge3J9CiMgMmEKQV9MU1IgPSBzb2x2ZSh0KGZ1bGxfVEMpICUqJSBmdWxsX1RDKSAlKiUgdChmdWxsX1RDKSAlKiUgc2NhbGVkX1gKRF9MU1IgPSBzY2FsZWRfWCUqJSBhYnModChBX0xTUikpCgoKIyAgUGxvdCBhbmQgaGVhdG1hcCB0b2dldGhlcgpwbG90KERfTFNSWzE6bnJvdyhEX0xTUiksMV0sIHR5cGUgPSAibGluZSIpCgojIHByb2R1Y2UgbGVnZW5kIGZvciBoZWF0bWFwCiNoZWF0bWFwKG1hdHJpeChBWzEsMTpuY29sKEFfTFNSKV0sIDIxLCAyMSksIG5jb2wsIENvbHYgPSBOQSwgUm93diA9IE5BLCBzY2FsZT0iY29sdW1uIikKCgoKZm9yIChudW0gaW4gMTo2KSB7CiAgZGYgPSBhcy5kYXRhLmZyYW1lKG1hdHJpeChyZXAoMCwgNDQxICogMyksIDQ0MSwgMykpCiAgY291bnQgPSAxCiAgY29sID0gYWJzKEFfTFNSW251bSxdKQogIGZvciAoaSBpbiAxOjIxKSB7CiAgICBmb3IgKGogaW4gMToyMSkgewogICAgICBkZltjb3VudCwgMV0gPSBpCiAgICAgIGRmW2NvdW50LCAyXSA9IGoKICAgICAgZGZbY291bnQsIDNdID0gY29sWyhpLTEpICogMjEgKyBqXQogICAgICBjb3VudCA9IGNvdW50ICsgMQogICAgfQogICAgCiAgfQogIGdncGxvdChkZiwgYWVzKFYxLCBWMiwgZmlsbD0gVjMpKSArIAogIGdlb21fdGlsZSgpICsgc2NhbGVfZmlsbF9ncmFkaWVudChsb3c9ImJsdWUiLCBoaWdoPSJyZWQiKSArCiAgdGhlbWVfaXBzdW0oKQp9CgoKCgoKCiMgV0h5IGxpbmVhciByZWxvIGJldCAzcmQgY29sIGRfbHNyIGFuZCAzMHRoIG9mIFgKIyBTZWNvbmRseSBpZiB5b3UgbG9vayBhdCB0aGUgc2xpY2UgbnVtYmVyIDMgYnVpbGQgdXAgYnkgdGhpcmQgU00uIE5vdyByZWFkIHRoaXMgM3JkIFNNIGNvbHVtbi13aXNlLCB5b3Ugd2lsbCBmaW5kIHRoYXQgMzB0aCAoMm5kIGNvbHVtbiA5dGggcm93ID0zMHRoIHBpeGVsKSBwaXhlbCBwb3NpdGlvbiB3YXMgZmlsbGVkIGJ5IHRoaXMgM3JkIFNNIGFuZCB0aGVyZWZvcmUgdGhpcmQgVEMgaXMgdGhlIG9ubHkgdGltZSBjb3Vyc2UgdGhhdCBjb25zdHJ1Y3RzIDMwdGggY29sdW1uIG9mIFggYW5kIHRoaXMgaXMgdGhlIFRDIHRoYXQgeW91IHJldHJpZXZlIGluIDNyZCBEbHNyLiBUaGVyZSBtaWdodCBiZSBzb21lIHZhcmlhbmNlIGluIHRoZSBzY2F0dGVyIHBsb3QgYnV0IHlvdSB3aWxsIGJlIGFibGUgdG8gY2FwdHVyZSB2aXN1YWxseSB0aGUgbGluZWFyIHRyZW5kIGJldHdlZW4gM3JkIERsc3IgYW5kIDMwdGggWC4KI2xpbmVhciByZWxvIC0gCnBsb3QoRF9MU1JbMTpucm93KERfTFNSKSwgM10sIHNjYWxlZF9YWzE6bnJvdyhzY2FsZWRfWCksIDMwXSwgeGxhYiA9ICIzcmQgY29sdW1uIG9mIEQoTFNSKSIsIHlsYWIgPSAiMzB0aCBjb2x1bW4gb2Ygc2NhbGVkIFgiKQoKIyBObyBsaW5lYXIgcmVsbwpwbG90KERfTFNSWzE6bnJvdyhEX0xTUiksIDRdLCBzY2FsZWRfWFsxOm5yb3coc2NhbGVkX1gpLCAzMF0sIHhsYWIgPSAiM3JkIGNvbHVtbiBvZiBEKExTUikiLCB5bGFiID0gIjMwdGggY29sdW1uIG9mIHNjYWxlZCBYIikKYGBgCmBgYHtyfQojIEdldHMgdGhlIHN1bSBvZiB0aGUgbWF4aW11bSBhYnNvbHV0ZSB2YWx1ZXMgaW4gdGhlIGNvcnJlbGF0aW9uIG1hdHJpeAptYXhfYWJzX2NvcnJfc3VtID0gZnVuY3Rpb24obWF0MSwgbWF0MikgewogIGNvcnJlX21hdCA9IGNvcihtYXQxLCBtYXQyKQogIGMgPSBnZXRfbWF4X2NvcnJlbGF0aW9ucyhjb3JyZV9tYXQpCiAgcmV0dXJuKHN1bShjKSkKfQoKIyBSZXR1cm5zIHRoZSB2YWx1ZXMgb2YgdGhlIG1heGltdW0gYWJzb2x1dGUgdmFsdWVzIGZvciBlYWNoIHJvdwpnZXRfbWF4X2NvcnJlbGF0aW9ucyA9IGZ1bmN0aW9uKGNvcl9tYXQpIHsKICBjX3ZhciA9IGMocmVwKDAsNikpCiAgYWJzX21hdCA9IGFicyhjb3JfbWF0KQogIGNvbHMgPSBtYXguY29sKGFic19tYXQpCiAgCiAgZm9yIChpIGluIDE6NikgewogICAgY192YXJbaV0gPSBjb3JfbWF0W2ksIGNvbHNbaV1dCiAgfQogIHJldHVybihjX3ZhcikKfQpgYGAKCgpgYGB7cn0KIzJiClYgPSA0NDEKTVNFID0gbWF0cml4KDAsIDEsIDIpCmZvciAobnVtIGluIHNlcShmcm9tID0gMCwgdG8gPSAxLCBieSA9IDAuMDEpKSB7CiAgbGFtZGEgPSBudW0KICBsYW1kYV9oYXQgPSBsYW1kYSAqIFYKICBBX1JSID0gc29sdmUodChmdWxsX1RDKSAlKiUgZnVsbF9UQyArIGRpYWcobnJvdyh0KGZ1bGxfVEMpKSkgKiBsYW1kYV9oYXQpICUqJSB0KGZ1bGxfVEMpICUqJSBzY2FsZWRfWAogIG1zZSA9IChzdW0oKHNjYWxlZF9YIC0gZnVsbF9UQyAlKiUgYWJzKEFfUlIpKV4yKSkgKyBsYW1kYV9oYXQgKiBzdW0oQV9SUl4yKQogIE1TRSA9IHJiaW5kKE1TRSwgYyhsYW1kYSwgbXNlKSkKfQoKcGxvdChNU0VbMjpucm93KE1TRSksIDFdLCBNU0VbMjpucm93KE1TRSksIDJdKQoKIyBGcm9tIHRoZSBncmFwaCBpdCBsb29rcyBsaWtlIGJlc3QgdmFsdWUgb2YgbGFtZGEgaXMgMQpBX1JSID0gc29sdmUodChmdWxsX1RDKSAlKiUgZnVsbF9UQyArIGRpYWcobnJvdyh0KGZ1bGxfVEMpKSkgKiAxKiBWKSAlKiUgdChmdWxsX1RDKSAlKiUgc2NhbGVkX1gKRF9SUiA9IHNjYWxlZF9YJSolIGFicyh0KEFfUlIpKQpjX3Rsc3IgPSBtYXhfYWJzX2NvcnJfc3VtKGZ1bGxfVEMsIERfTFNSKQpjX3RyciA9IG1heF9hYnNfY29ycl9zdW0oZnVsbF9UQywgRF9SUikKCiMgNC40MDcyNDEKY190bHNyCgojIDQuNDc2NTMKY190cnIKCgpjX3RyciAtIGNfdGxzcgoKCmBgYApgYGB7cn0KCmxhbWRhID0gMTAwMApsYW1kYV9oYXQgPSBsYW1kYSAqIFYKQV8xMDAwID0gc29sdmUodChmdWxsX1RDKSAlKiUgZnVsbF9UQyArIGRpYWcobnJvdyh0KGZ1bGxfVEMpKSkgKiBsYW1kYV9oYXQpICUqJSB0KGZ1bGxfVEMpICUqJSBzY2FsZWRfWAoKCiMgQXMgQV9MU1IgZGVjcmVhc2UsIEFfUlIgdmFsdWVzIGRlY3JlYXNlCnBsb3QoYWJzKEFfMTAwMFsxOm5yb3coQV8xMDAwKSwgMV0pLCBhYnMoQV9MU1JbMTpucm93KEFfTFNSKSwgMV0pLCB4bGFiID0gIkFfUlIiLCB5bGFiID0gIkFfTFNSIikKCmBgYAoKYGBge3J9CmNyZWF0ZV9ub2lzZV9tYXRyaXggPSBmdW5jdGlvbihzZWVkKSB7CiAgc2V0LnNlZWQoc2VlZCA9IHNlZWQpCiAgbm9pc2VfdCA9IHJub3JtKDYgKiAyNDAsIG1lYW4gPSAwLCBzZCA9IHNxcnQoMC4yNSkpCiAgbm9pc2VfcyA9IHJub3JtKDYgKiA0NDEsIG1lYW4gPSAwLCBzZCA9IHNxcnQoMC4xNSkpCgoKICBub2lzZV9zID0gbWF0cml4KG5vaXNlX3MsIDYsIDQ0MSkKICBub2lzZV90ID0gbWF0cml4KG5vaXNlX3QsIDI0MCwgNikKICAKICBYID0gKGZ1bGxfVEMgKyBub2lzZV90KSAlKiUgKHRtcFNNICsgbm9pc2VfcykKICBuZXdfWCA9IHNjYWxlKFgpCiAgcmV0dXJuKG5ld19YKQp9Cm5vaXNlXzEgPSBjcmVhdGVfbm9pc2VfbWF0cml4KDMwMDM1KQpub2lzZV8yID0gY3JlYXRlX25vaXNlX21hdHJpeCgzMDAzNykKbm9pc2VfMyA9IGNyZWF0ZV9ub2lzZV9tYXRyaXgoMzAwMzgpCm5vaXNlXzQgPSBjcmVhdGVfbm9pc2VfbWF0cml4KDMwMDM2KQpub2lzZV81ID0gY3JlYXRlX25vaXNlX21hdHJpeCgzMDAzOSkKbm9pc2VfNiA9IGNyZWF0ZV9ub2lzZV9tYXRyaXgoMzAwNDApCm5vaXNlXzcgPSBjcmVhdGVfbm9pc2VfbWF0cml4KDMwMDQxKQpub2lzZV84ID0gY3JlYXRlX25vaXNlX21hdHJpeCgzMDA0MikKbm9pc2VfOSA9IGNyZWF0ZV9ub2lzZV9tYXRyaXgoMzAwNDMpCm5vaXNlXzEwID0gY3JlYXRlX25vaXNlX21hdHJpeCgzMDA0NCkKCgpub2lzZV9kZiA9IGxpc3Qobm9pc2VfMSwgbm9pc2VfOSwgbm9pc2VfOCwgbm9pc2VfNywgbm9pc2VfNiwgbm9pc2VfNSwKICAgICAgICAgICAgICAgIG5vaXNlXzIsIG5vaXNlXzMsIG5vaXNlXzQsIG5vaXNlXzEwCiAgICAgICAgICAgICAgICApCiNmb3IobWF0IGluIG5vaXNlX2RmKSB7CiMgIHByaW50KG1hdFsxOjUsIDE6M10KI30KCgpgYGAKCmBgYHtyfQojMmMKTFJfTVNFID0gbWF0cml4KE5BLCAyMSwgMTApCk4gPSAyNDAKeDEgPSAyMQp4MiA9IDIxCm5zcmNzID0gNgpWID0geDEgKiB4MgppID0gMQpqID0gMQpzdGVwIDwtIDEvKG5vcm0oZnVsbF9UQyAlKiUgdChmdWxsX1RDKSkgKiAxLjEpCgpmb3IgKG51bSBpbiBzZXEoMCwxLDAuMDUpKSB7CiAgZm9yICh4X3JlcCBpbiBub2lzZV9kZikgewogICAgdGhyIDwtIG51bSpOKnN0ZXAKICAgIEFvIDwtIG1hdHJpeCgwLCBuc3JjcywgMSkKICAgIEEgPC0gbWF0cml4KDAsIG5zcmNzLCAxKQogICAgQWxyIDwtIG1hdHJpeCgwLCBuc3JjcywgeDEqeDIpCiAgICBmb3IgKGsgaW4gMTooeDEqeDIpKSB7CiAgICAgIEEgPC0gQW8rc3RlcCoodChmdWxsX1RDKSAlKiUgKHhfcmVwWyxrXS0oZnVsbF9UQyUqJUFvKSkpCiAgICAgIEEgPC0gKDEvKDErdGhyKSkgKiAoc2lnbihBKSpwbWF4KHJlcGxpY2F0ZShuc3JjcywgMCksIGFicyhBKS10aHIpKQogICAgICAKICAgICAgZm9yIChyZXAgaW4gMToxMCkgewogICAgICAgIEFvIDwtIEEKICAgICAgICBBIDwtIEFvK3N0ZXAgKiAodChmdWxsX1RDKSUqJSh4X3JlcFssa10tKGZ1bGxfVEMlKiVBbykpKQogICAgICAgIEEgPC0gKDEvKDErdGhyKSkgKiAoc2lnbihBKSpwbWF4KHJlcGxpY2F0ZShuc3JjcywgMCksIGFicyhBKS10aHIpKQogICAgICB9CiAgICAgIAogICAgICBBbHJbLGtdIDwtIEEKICAgIH0KICAgIERsciA9IHhfcmVwICUqJSBhYnModChBbHIpKQogICAgTFJfTVNFW2ksal0gPSBzdW0oc3VtKCh4X3JlcC1EbHIlKiVhYnMoQWxyKSleMikpIC8gKE4qVikKICAgIGogPSBqICsgMQogIH0KICBqID0gMQogIGkgPSBpICsgMQp9CgpMUl9hdmdzID0gcm93TWVhbnMoTFJfTVNFKQpwbG90KCBzZXEoMCwxLCAwLjA1KSwgTFJfYXZncywgeWxhYiA9ICJNU0UoTFIpIikKIyBNaW4gdmFsdWUgaXMgMC41ODAgd2hlbiBwID0gMC41NSwgYXQgc2FtZSB2YWx1ZSBpdCBkaXZlcmdlcwoKIyBTaW5jZSBmaXJzdCBpbmRleCBpcyBhdCAwLCB3ZSB3YW50IHRoZSAxMnRoIHZhbHVlIGJldHdlZW4gMCBhbmQgMSB3aGljaCBpcyAwLjYgCiMgU28gd2UgbmVlZCB0byBpbmRleCB0aGUgMTN0aCB2YWx1ZQppbmRleCA9IHdoaWNoLm1pbihMUl9hdmdzKQpyaG8gPSAoaW5kZXggLSAxKSAqIDAuMDUKYGBgCgpgYGB7cn0KIyAyZAoKdGhyIDwtIHJobypOKnN0ZXAKQW8gPC0gbWF0cml4KDAsIG5zcmNzLCAxKQpBIDwtIG1hdHJpeCgwLCBuc3JjcywgMSkKQWxyIDwtIG1hdHJpeCgwLCBuc3JjcywgeDEqeDIpCmZvciAoayBpbiAxOih4MSp4MikpIHsKICBBIDwtIEFvK3N0ZXAqKHQoZnVsbF9UQykgJSolIChzY2FsZWRfWFssa10tKGZ1bGxfVEMlKiVBbykpKQogIEEgPC0gKDEvKDErdGhyKSkgKiAoc2lnbihBKSpwbWF4KHJlcGxpY2F0ZShuc3JjcywgMCksIGFicyhBKS10aHIpKQogIAogIGZvciAocmVwIGluIDE6MTApIHsKICAgIEFvIDwtIEEKICAgIEEgPC0gQW8rc3RlcCAqICh0KGZ1bGxfVEMpJSolKHNjYWxlZF9YWyxrXS0oZnVsbF9UQyUqJUFvKSkpCiAgICBBIDwtICgxLygxK3RocikpICogKHNpZ24oQSkqcG1heChyZXBsaWNhdGUobnNyY3MsIDApLCBhYnMoQSktdGhyKSkKICB9CiAgCiAgQWxyWyxrXSA8LSBBCn0KRGxyID0gc2NhbGVkX1ggJSolIGFicyh0KEFscikpCmBgYAoKYGBge3J9CmN0cnIgPSBtYXhfYWJzX2NvcnJfc3VtKGZ1bGxfVEMsIERfUlIpCmNzcnIgPSBtYXhfYWJzX2NvcnJfc3VtKHQodG1wU00pLCB0KEFfUlIpKQpjdGxyID0gbWF4X2Fic19jb3JyX3N1bShmdWxsX1RDLCBEbHIpCmNzbHI9IG1heF9hYnNfY29ycl9zdW0odCh0bXBTTSksIHQoQWxyKSkKCiMgQm90aCBhcmUgdHJ1ZQpjdGxyIC0gY3RyciAjPSAwLjc0MTM0ODcKY3NsciAtIGNzcnIgIz0gMC42NTk1MTUxCmN0bHIgIz0gNS4yMTc4NzkKY3RyciAjPSA0LjQ3NjUzCmNzbHIgIz0gMi45NTMxMDEKY3NyciAjPSAyLjI5MzU4NgoKcGxvdChEX1JSWywgMTo0XSkKcGxvdChEbHJbLCAxOjRdKQpwbG90KEFsclsxOjQsXSkKcGxvdChBX1JSWzE6NCwgXSkKCnBhcihtZnJvdyA9IGMoMSwyKSkKZm9yIChudW0gaW4gMTo2KSB7CiAgcGxvdChEX1JSWywgNl0sIHR5cGUgPSAibGluZSIpCiAgcGxvdChEbHJbLCA2XSwgdHlwZSA9ICJsaW5lIikKfQoKCgpkZiA9IGFzLmRhdGEuZnJhbWUobWF0cml4KHJlcCgwLCA0NDEgKiAzKSwgNDQxLCAzKSkKY291bnQgPSAxCnZlYyA9IGFicyhBX1JSW251bSxdKQpmb3IgKG51bSBpbiAxOjYpIHsKICBmb3IgKGkgaW4gMToyMSkgewogICAgZm9yIChqIGluIDE6MjEpIHsKICAgICAgZGZbY291bnQsIDFdID0gaQogICAgICBkZltjb3VudCwgMl0gPSBqCiAgICAgIGRmW2NvdW50LCAzXSA9IHZlY1soaS0xKSAqIDIxICsgal0KICAgICAgY291bnQgPSBjb3VudCArIDEKICAgIH0KICAKICB9Cn0KCmRmID0gYXMuZGF0YS5mcmFtZShtYXRyaXgocmVwKDAsIDQ0MSAqIDMpLCA0NDEsIDMpKQpjb3VudCA9IDEKdmVjID0gYWJzKEFscltudW0sXSkKZm9yIChudW0gaW4gMTo2KSB7CiAgZm9yIChpIGluIDE6MjEpIHsKICAgIGZvciAoaiBpbiAxOjIxKSB7CiAgICAgIGRmW2NvdW50LCAxXSA9IGkKICAgICAgZGZbY291bnQsIDJdID0gagogICAgICBkZltjb3VudCwgM10gPSB2ZWNbKGktMSkgKiAyMSArIGpdCiAgICAgIGNvdW50ID0gY291bnQgKyAxCiAgICB9CiAgCiAgfQp9CgpnZ3Bsb3QoZGYsIGFlcyhWMSwgVjIsIGZpbGw9IFYzKSkgKyAKICBnZW9tX3RpbGUoKSArIHNjYWxlX2ZpbGxfZ3JhZGllbnQobG93PSJibHVlIiwgaGlnaD0icmVkIikgKwogIHRoZW1lX2lwc3VtKCkKCmBgYAoKYGBge3J9CnBhcihtZnJvdz1jKDEsMikpCiMgMmUKI3N2ZChmdWxsX1RDKQpwbG90KHN2ZChmdWxsX1RDKSRkLCB5bGFiID0gIkVpZ2VudmFsdWUiKQojUEMgNiBzbWFsbGVzdAoKcGxvdChzdmQoZnVsbF9UQykkdVsxOjI0MCwgNl0pCnBsb3QoZnVsbF9UQ1sxOjI0MCwgNl0pCgpwbG90KHN2ZChmdWxsX1RDKSR1WzE6MjQwLCA1XSkKcGxvdChmdWxsX1RDWzE6MjQwLCA1XSkKCnBsb3Qoc3ZkKGZ1bGxfVEMpJHVbMToyNDAsIDRdKQpwbG90KGZ1bGxfVENbMToyNDAsIDRdKQoKcGxvdChzdmQoZnVsbF9UQykkdVsxOjI0MCwgM10pCnBsb3QoZnVsbF9UQ1sxOjI0MCwgM10pCgpwbG90KHN2ZChmdWxsX1RDKSR1WzE6MjQwLCAyXSkKcGxvdChmdWxsX1RDWzE6MjQwLCAyXSkKCnBsb3Qoc3ZkKGZ1bGxfVEMpJHVbMToyNDAsIDFdKQpwbG90KGZ1bGxfVENbMToyNDAsIDFdKQoKCmBgYAoKCgoK